
# set working directory
setwd("~/Research/Mycorrhizae/Fieldwork/Isotopes")

library(lme4)
library(nlme)
library(lmerTest)
library(MuMIn)
library(htmlTable)
library(ggplot2)
library(Rmisc)
library(SciViews)
library(dplyr)
library(MASS)
library(gridExtra)
library(ggeasy)
library(grid)
library(regclass)

# set scientific notation off
options(scipen=999)

# define predictive R^2 functions for later
PRESS <- function(linear.model) {
  #' calculate the predictive residuals
  pr <- residuals(linear.model)/(1-lm.influence(linear.model)$hat)
  #' calculate the PRESS
  PRESS <- sum(pr^2)
  
  return(PRESS)
}

pred_r_squared <- function(linear.model) {
  #' Use anova() to get the sum of squares for the linear model
  lm.anova <- anova(linear.model)
  #' Calculate the total sum of squares
  tss <- sum(lm.anova$'Sum Sq')
  # Calculate the predictive R^2
  pred.r.squared <- 1-PRESS(linear.model)/(tss)
  
  return(pred.r.squared)
}

#########################################
# CREATING GRAPHS OF RAW DATA AND MEANS
################################################

# upload data
fungi.data <- read.csv("WSU mushroom isotopes for analysis.csv")

# make distance numeric
fungi.data$Distance <- as.numeric(fungi.data$Distance)

# put distance 4 into 6 bin for analysis
for (i in 1:length(fungi.data$Distance)){
  if (fungi.data$Distance[i]==4){
    fungi.data$Distance[i] <- 6
  }
}

# by stream 
fungi.data <- fungi.data[ which(fungi.data$Stream=="Hansen"),]
fungi.data <- fungi.data[ which(fungi.data$Stream=="Happy"),]
fungi.data <- fungi.data[ which(fungi.data$Stream=="Yako"),]

# remove saprotrophs
fungi.data <- fungi.data[which(fungi.data$trophic.status != "saprotroph"),]
unique(fungi.data$Genus)

# by exploration type
fungi.data <- fungi.data[which(fungi.data$exploration.type == "short distance"),]
                         
# by hansen side
fungi.data <- fungi.data[ which(fungi.data$Side=="SL"),]
fungi.data <- fungi.data[ which(fungi.data$Side=="SR"),]

# by distance
fungi.data <- fungi.data[which(fungi.data$Distance < 21),]

# insert Russula genus guesses if needed
fungi.data$Genus[71] <- "Russula"
fungi.data$Genus[73] <- "Russula"
fungi.data$Genus[76] <- "Russula"
fungi.data$Genus[81] <- "Russula"
fungi.data$Genus[103] <- "Russula"
fungi.data$Genus[109] <- "Russula"

# by trophic status
fungi.data <- fungi.data[which(fungi.data$trophic.status == "mycorrhizal"),]

# remove those points with decomposing salmon carcass
fungi.data <- fungi.data[which(fungi.data$decomposing.carcass == "no"),]

# fungal species
fungi.data <- fungi.data[ which(fungi.data$Genus=="Lactarius"),]
fungi.data <- fungi.data[ which(fungi.data$Genus=="Russula"),]
fungi.data <- fungi.data[ which(fungi.data$Genus=="Cortinarius"),]
fungi.data <- fungi.data[ which(fungi.data$Genus=="Paxillus"),]
fungi.data <- fungi.data[ which(fungi.data$Genus=="Naucoria"),]

# choose transect
fungi.data <- fungi.data[ which(fungi.data$Transect=="S1"),]

# for total N graphs on Hansen left bank, remove the outlier Lactarius 
# at 100 m that has very high C:N
fungi.data <- fungi.data[which(fungi.data$total.N < 0.1),]

# examine by Hansen bank 
tgc <- summarySE(fungi.data, measurevar="d15N", groupvars=c("Distance", "Side"))
tgc

write.csv(file="fungi_mean_d15N.csv", x=tgc)

# for Happy remove averages for single data points
tgc <- tgc[-c(3,8),]
tgc

# for Yako left bank remove averages for single data points
tgc <- tgc[-c(1,7),]
tgc

# for Yako right bank remove averages for single data points
tgc <- tgc[-c(2,5),]
tgc

# for Hansen left bank remove averages for single data points
tgc <- tgc[-c(2,8:10),]

# for Hansen right bank remove averages for single data points
tgc <- tgc[-7,]

p1 <- ggplot(tgc[which(tgc$Side=="SL"),], aes(x=Distance, y=d15N)) + geom_point(position = position_dodge(.9), 
   stat = "summary", fun = "mean") + 
  geom_errorbar(aes(ymin=d15N-se, ymax=d15N+se), width=.2,position=position_dodge(.9)) +
  theme_classic() + 
  geom_point(data=fungi.data[which(fungi.data$Side=="SL"),], aes(x=Distance, y=d15N), 
             alpha=0.2) + ylim(-5,19) + scale_x_reverse(breaks=c(1,6,10,20, 40, 60, 100)) 

p2 <- ggplot(tgc[which(tgc$Side=="SR"),], aes(x=Distance, y=d15N)) + geom_point(position = position_dodge(.9), 
   stat = "summary", fun = "mean") + 
  geom_errorbar(aes(ymin=d15N-se, ymax=d15N+se), width=.2,position=position_dodge(.9)) +
  theme_classic() + 
  geom_point(data=fungi.data[which(fungi.data$Side=="SR"),], aes(x=Distance, y=d15N), 
             alpha=0.2) + ylim(-5,19) + scale_x_continuous(breaks=c(1,6,10,20, 40, 60, 100)) 

grid.arrange(p1, p2, nrow=1)

ggsave(plot=p, width=5, height=3, dpi=300, filename = "d15N_raw_fungi_Hansen_SR1_sample.jpg")

# mean and errors
p <- ggplot(fungi.data, aes(x=Side, y=d15N)) + 
  geom_errorbar(aes(ymin=d15N-se, ymax=d15N+se), width=.1) + geom_line() + 
  geom_point() + theme_bw() + scale_x_discrete(labels=c('Enhanced', 'Depleted')) +
  ylab("Fungal sporocarp d15N") 
p

t.test(fungi.data$d15N[which(fungi.data$Side=="SL")], 
       fungi.data$d15N[which(fungi.data$Side=="SR")])

ggsave(plot=p, width=3, height=3, dpi=300, filename = "dC13_Hansen_fungi_byside.jpg")

# boxplot
ggplot(fungi.data, aes(x=Side, y=d15N, color=Side)) + geom_boxplot() + 
  stat_summary(fun.y=mean, geom="point", shape=8, size=4) + 
  scale_color_brewer(palette="Dark2") +
  theme_classic() +
  theme(axis.title.x = element_text(size = 18),
        axis.title.y = element_text(size = 18),
        axis.text.x = element_text(size = 14),
        axis.text.y = element_text(size = 14), legend.position="none")

# Hansen
# for C:N, use ylim(4,14)
# for d15N, use ylim(-5,19)
# for %N, use ylim(2.5, 10)
# for d13C, use ylim(-30,-21)
# for total N, use ylim(0,0.044), UNLESS want to include outlier then ylim(0,0.132)

# for Lactarius dN15, use ylim(2.5,15)
# for Russula dN15, use ylim(0,12)
# for Cortinarius dN15, use ylim (5.5, 17)

# Yako
# for d15N, use ylim(0.025, 20)

# examine at Happy
tgc <- summarySE(fungi.data, measurevar="d13C", groupvars=c("Distance","Species"))
ggplot(tgc, aes(x=Distance, y=d13C, color=Species)) + geom_point(position = position_dodge(.9), 
                                                                 stat = "summary", fun = "mean") + geom_line() +
  geom_errorbar(aes(ymin=d13C-se, ymax=d13C+se), width=.2,position=position_dodge(.9)) +
  theme_classic() + 
  geom_point(data=fungi.data, aes(x=Distance, y=d13C), alpha=0.2)  +
  ylim(-34,-29)

# Happy
# for dC13, use ylim(-33,-28)
# for C:N, use ylim(15,56)
# Yako
# for C:N, use ylim(17,55)
# for d15N, use ylim(-4,6)

# C:N graph
fungi.data <- read.csv("WSU mushroom isotopes for analysis.csv")

# if using full dataset, remove Aleknagik data
fungi.data <- fungi.data[which(fungi.data$Stream != "Aleknagik"),]
p <- ggplot(fungi.data, aes(x=X.C, y=X.N, color=trophic.status)) + geom_point(size=3) + theme_classic() + 
  theme(axis.title.x = element_text(size = 30),
        axis.title.y = element_text(size = 30),
        axis.text.x = element_text(size = 26),
        axis.text.y = element_text(size = 26), plot.title = element_text(hjust = 0.5),
        legend.title = element_blank(),
        legend.text = element_text(size = 26)) +
  xlab("%C") + ylab("%N") 

ggsave(plot=p, width=15, height=9, dpi=300, filename = "Fig4.jpg")

# dC13 vs d15N graph
fungi.data <- read.csv("WSU mushroom isotopes for analysis.csv")

# if using full dataset, remove Aleknagik data
fungi.data <- fungi.data[which(fungi.data$Stream != "Aleknagik"),]
p <- ggplot(fungi.data, aes(x=dC13, y=d15N, color=trophic.status)) + geom_point(size=3) + theme_classic() + 
  theme(axis.title.x = element_text(size = 30),
        axis.title.y = element_text(size = 30),
        axis.text.x = element_text(size = 26),
        axis.text.y = element_text(size = 26), plot.title = element_text(hjust = 0.5),
        legend.title = element_blank(),
        legend.text = element_text(size = 26)) + 
  labs(x=expression(delta^13*"C (\u2030 vs. VPDB)"), y=expression(delta^15*"N (\u2030 vs. air)"))

p   

ggsave(plot=p, width=15, height=9, dpi=300, filename = "Fig5.jpg")


#######################################################################
# DATA PREP FOR MODEL ANALYSIS BELOW 
#####################################################################
fungi.data <- read.csv("WSU mushroom isotopes for analysis.csv")

# make distance numeric
fungi.data$Distance <- as.numeric(fungi.data$Distance)

# by stream 
fungi.data <- fungi.data[ which(fungi.data$Stream=="Hansen"),]

# select only 1-20 m  
fungi.data <- fungi.data[ which(fungi.data$Distance < 7),]

#######################################################################
# C:N, dN15 and dC13 by distance and bank at Hansen
#####################################################################

##############################################
# C:N #
#########################################
# select data
fungi.data <- read.csv("WSU mushroom isotopes for analysis.csv")

# make distance numeric
fungi.data$Distance <- as.numeric(fungi.data$Distance)

# choose trophic status (unknown species are likely saprotrophs)
fungi.data <- fungi.data[which(fungi.data$trophic.status == "mycorrhizal" | fungi.data$trophic.status=="unknown"),]
fungi.data <- fungi.data[which(fungi.data$trophic.status == "saprotroph"),]

# choose exploration type
fungi.data <- fungi.data[which(fungi.data$exploration.type == "long distance"),]
fungi.data <- fungi.data[which(fungi.data$exploration.type == "medium distance smooth"),]
fungi.data <- fungi.data[which(fungi.data$exploration.type == "medium distance fringe"),]
fungi.data <- fungi.data[which(fungi.data$exploration.type == "short distance"),]

# by stream 
#fungi.data <- fungi.data[ which(fungi.data$Stream=="Hansen"),]

# select only 1-20 m  
#fungi.data <- fungi.data[ which(fungi.data$Distance < 21),]

# scale variables
test1 <- scale(fungi.data$organic.GWC, scale=FALSE)
fungi.data$organic.GWC <- test1[,1]

test2 <- scale(fungi.data$mineral.GWC, scale=FALSE)
fungi.data$mineral.GWC <- test2[,1]

test3 <- scale(fungi.data$organic.soil.C.N, scale=FALSE)
fungi.data$organic.soil.C.N <- test3[,1]

test4 <- scale(fungi.data$mineral.soil.C.N, scale=FALSE)
fungi.data$mineral.soil.C.N <- test4[,1]

fungi.data$decomposing.carcass <- as.factor(fungi.data$decomposing.carcass)
fungi.data$carcass.deposition <- as.factor(fungi.data$carcass.deposition)
fungi.data$trophic.status <- as.factor(fungi.data$trophic.status)
fungi.data$carcass.deposition.ord <- as.factor(fungi.data$carcass.deposition.ord)
fungi.data$carcass.load <- as.numeric(fungi.data$carcass.load)

# if using full dataset, remove Aleknagik data
fungi.data <- fungi.data[which(fungi.data$Stream != "Aleknagik"),]

# remove NA values for organic and mineral soil dN15 and C:N
library(tidyr)
fungi.data <- fungi.data[!is.na(fungi.data$organic.soil.C.N),] 
fungi.data <- fungi.data[!is.na(fungi.data$mineral.soil.C.N),] 

# if using foilar N:P, remove any rows that do not have value for foliar N:P
#fungi.data <- fungi.data[!is.na(fungi.data$paper.birch.N.P),]

# remove certain variables (lose 13 data points for fungi)
#fungi.data <- fungi.data[which(fungi.data$trophic.status == "mycorrhizal"),]

# if including plant C:N, need to drop a few rows that did not have data for plant C:N 
fungi.data <- fungi.data %>% drop_na(white.spruce.C.N)
fungi.data <- fungi.data %>% drop_na(paper.birch.C.N)

#fungi.data <- fungi.data[which(fungi.data$Genus=="Laccaria"),]

# parameters for BOTH SIDES
options(na.action = "na.fail")
all.parms<-lm(log(C.N) ~ 
                #Side + 
                ln(Distance) + 
                #Side:ln(Distance) +
                #Side:I(ln(Distance)^2) + 
                #Side:I(ln(Distance)^3) + 
                #Side:I(ln(Distance)^4) + 
                paper.birch.C.N + 
                white.spruce.C.N +
                organic.GWC + 
                mineral.GWC +
                organic.soil.C.N + 
                mineral.soil.C.N +
                carcass.load + 
                #white.spruce.N.P + 
                #paper.birch.N.P +
                trophic.status + 
                carcass.deposition +
                #carcass.deposition.ord + 
                decomposing.carcass,
              data=fungi.data)
results<-dredge(all.parms)
subset(results, delta <2)
subset(results, delta == 0)

dredge.model <- lmer(formula = log(C.N) ~ carcass.deposition + carcass.load + paper.birch.C.N  + white.spruce.C.N +
                       trophic.status + (1|Genus) + (1|Stream:Transect), data=fungi.data)
summary(dredge.model)

# test dredge model
dredge.model <- lm(formula = log(C.N) ~ carcass.deposition + carcass.load + paper.birch.C.N  + white.spruce.C.N +
                       trophic.status, data=fungi.data)
summary(dredge.model)

# step model
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

# step model with random effect
step.model <- lm(log(C.N) ~ organic.GWC + organic.soil.C.N + carcass.deposition, data = fungi.data)
summary(step.model)

# for only mycorrhizal fungi
dredge.model <- lmer(C.N ~ organic.GWC + ln(Distance) + mineral.GWC + organic.soil.C.N +
                     paper.birch.C.N + white.spruce.C.N + (1|Genus) + 
                       (1|Stream:Transect), data=fungi.data)
summary(dredge.model)

# all data
dredge.model <- lmer(C.N ~ ln(Distance) + paper.birch.C.N + white.spruce.C.N + organic.GWC + mineral.GWC + organic.soil.C.N + 
                     mineral.soil.C.N + trophic.status + carcass.deposition + decomposing.carcass + carcass.load + (1|Genus) + 
                     (1|Stream:Transect), data=fungi.data)
summary(dredge.model)

# model for single genera
dredge.model <- lm(C.N ~ carcass.deposition + decomposing.carcass, data=fungi.data)
summary(dredge.model)

dredge.model <- lm(C.N ~ carcass.deposition.ord + paper.birch.C.N + trophic.status + white.spruce.C.N, data=fungi.data)
summary(dredge.model)

# model for all generas from 1-100
dredge.model <- lm(C.N ~ carcass.deposition.ord + ln(Distance) + paper.birch.C.N + trophic.status + white.spruce.C.N, data=fungi.data)
summary(dredge.model)

dredge.model <- lmer(C.N ~ carcass.deposition + ln(Distance) + paper.birch.C.N + trophic.status + white.spruce.C.N + 
                       (1|Genus), data=fungi.data)
summary(dredge.model)

# get conditional r^2 (which takes into account random effects, second number)
r.squaredGLMM(lme4::lmer(C.N ~ carcass.deposition + carcass.load + paper.birch.C.N +  
                                  white.spruce.C.N + trophic.status + (1 | Genus) + (1 | Stream:Transect), data=fungi.data))


# get percentage variance explained by each variable using anova
af <- anova(dredge.model)
afss <- af$"Sum Sq"
blah <- cbind(af,PctExp=afss/sum(afss)*100)
blah$PctExp

# calculate VIF
library(car)
car::vif(dredge.model)

# for fungal 1:20 SL
dredge.model <- lm(C.N ~ carcass.deposition + trophic.status + ln(Distance) + organic.GWC, data=fungi.data)
summary(dredge.model)

# for 1-20
dredge.model <- lm(C.N ~ Side + Side:I(ln(Distance)^3), data=fungi.data)
summary(dredge.model)

# step model
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

# for 1-100, both dredge and step models give very low predictive R squared
# looking at the raw data, it is pretty variable, especially on the left side

pred_r_squared(dredge.model)

# examine effect of distance on C:N ratio
mod1 <- lm(C.N ~ Side, data=fungi.data)
mod2 <- lm(C.N ~ Side + organic.GWC, data=fungi.data)
mod3 <- lm(C.N ~ Side + mineral.GWC, data=fungi.data)
mod4 <- lm(C.N ~ Side + organic.soil.C.N, data=fungi.data)
mod5 <- lm(C.N ~ Side + mineral.soil.C.N, data=fungi.data)
mod6 <- lm(C.N ~ ln(Distance) + organic.GWC, data=fungi.data)
mod7 <- lm(C.N ~ ln(Distance) + mineral.GWC, data=fungi.data)
mod8 <- lm(C.N ~ ln(Distance) + organic.soil.C.N, data=fungi.data)
mod9 <- lm(C.N ~ ln(Distance) + mineral.soil.C.N, data=fungi.data)
mod6a <- lm(C.N ~ I(ln(Distance)^2) + organic.GWC, data=fungi.data)
mod7a <- lm(C.N ~ I(ln(Distance)^2) + mineral.GWC, data=fungi.data)
mod8a <- lm(C.N ~ I(ln(Distance)^2) + organic.soil.C.N, data=fungi.data)
mod9a <- lm(C.N ~ I(ln(Distance)^2) + mineral.soil.C.N, data=fungi.data)
mod10a <- lm(C.N ~ ln(Distance), data=fungi.data)
mod10 <- lm(C.N ~ I(ln(Distance)^2), data=fungi.data)
mod11 <- lm(C.N ~ ln(Distance) + Side, data=fungi.data)
mod12 <- lm(C.N ~ ln(Distance) + Side + organic.GWC, data=fungi.data)
mod13 <- lm(C.N ~ ln(Distance) + Side + mineral.GWC, data=fungi.data)
mod14 <- lm(C.N ~ ln(Distance) + Side + organic.soil.C.N, data=fungi.data)
mod15 <- lm(C.N ~ ln(Distance) + Side + mineral.soil.C.N, data=fungi.data)
mod16 <- lm(C.N ~ organic.GWC, data=fungi.data)
mod16a <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance), data=fungi.data)
mod17 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + organic.GWC, data=fungi.data)
mod18 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + mineral.GWC, data=fungi.data)
mod19 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + organic.GWC + organic.soil.C.N, data=fungi.data)
mod20 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2), data=fungi.data)
mod21 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance), data=fungi.data)
mod22 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + organic.GWC, data=fungi.data)
mod23 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + mineral.GWC, data=fungi.data)
mod24 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2), data=fungi.data)
mod25 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + organic.GWC, data=fungi.data)
mod26 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + mineral.GWC, data=fungi.data)
mod27 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + organic.soil.C.N + organic.GWC + mineral.soil.C.N + mineral.GWC, data=fungi.data)
mod28 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3), data=fungi.data)
mod29 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.GWC, data=fungi.data)
mod30 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.GWC, data=fungi.data)
mod30a <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.soil.C.N, data=fungi.data)
mod31 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + mineral.GWC, data=fungi.data)
mod32 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.soil.C.N, data=fungi.data)
mod33 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.soil.C.N + organic.GWC, data=fungi.data)
mod34 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + mineral.soil.C.N + mineral.GWC, data=fungi.data)
mod35 <- lm(C.N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.soil.C.N + organic.GWC + mineral.soil.C.N + mineral.GWC, data=fungi.data)

# although model 15s have as low AIC as model 12, they do not include individual
# variables of side and distance which are important 
mod36 <- lm(C.N ~ Side:I(ln(Distance)^2), data=fungi.data)
mod37 <- lm(C.N ~ Side:I(ln(Distance)^3), data=fungi.data)
mod38 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE), data=fungi.data)
mod39 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE), data=fungi.data)
mod40 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC, data=fungi.data)
mod41 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC + organic.soil.C.N, data=fungi.data)
mod42 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE) + mineral.soil.C.N, data=fungi.data)
mod43 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC, data=fungi.data)
mod44 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.soil.C.N, data=fungi.data)
mod45 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC + organic.soil.C.N + mineral.soil.C.N + mineral.GWC, data=fungi.data)
mod46 <- lm(C.N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.soil.C.N + organic.GWC, data=fungi.data)
mod47 <- lm(C.N ~ Side:poly(ln(Distance), 4, raw=TRUE), data=fungi.data)
mod48 <- lm(C.N ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.GWC + organic.soil.C.N + mineral.soil.C.N + mineral.GWC, data=fungi.data)
mod49 <- lm(C.N ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.GWC + organic.soil.C.N, data=fungi.data)
mod50 <- lm(C.N ~ Side:poly(ln(Distance), 5, raw=TRUE), data=fungi.data)
mod51 <- lm(C.N ~ Side:poly(ln(Distance), 5, raw=TRUE) + organic.GWC + organic.soil.C.N + mineral.soil.C.N + mineral.GWC, data=fungi.data)
mod52 <- lm(C.N ~ Side:poly(ln(Distance), 5, raw=TRUE) + organic.GWC + organic.soil.C.N, data=fungi.data)

results <- AIC(mod1, mod2, mod3, mod4, mod5, mod6, mod7, mod8, mod9, mod10, 
               mod11, mod12, 
               mod13, mod14, mod15, mod16, mod16a, mod17, mod18, mod19, mod20, mod21, mod22, mod23, 
               mod24, mod25, mod26, mod27, mod28, mod29, mod30, mod30a, mod31, mod32, mod33, mod34,
               mod35, mod36, mod37, mod38, mod39, mod40, mod41, mod42, mod43, mod44, 
               mod45, mod46, mod47, mod48, mod49, mod50, mod51, mod52, step.model, dredge.model)
results <- AIC(mod6, mod7, mod8, mod9, mod10, mod16, mod6a, mod7a, mod8a, mod9a, mod10a, step.model, dredge.model, test.model)
ordered <- results[order(-results$AIC),]
ordered
summary(mod7)

# FUNGAL C:N 1-100
# some competing models, with dredge, step, mod7 and mod24. Dredge is best but no 
# difference between banks, step is a little overfitted with large humps on both 
# sides, and mod24 is opposing patterns close to banks. After examining close to 
# bank patterns, CHOSE MOD24 because same as for close bank and nicely fits pattern

# FUNGAL C:N 1-20
# top models include step and dredge

# select the model that will be used to plot predictions
select.mod <- dredge.model

# plot raw data with model estimates
pred_interval <- seq(1,20,0.2)
which.side <- "SL"
example <- data.frame(Distance=pred_interval, Side=which.side, 
          carcass.deposition = "yes",            
          total.N = mean(fungi.data$total.N[which(fungi.data$Side==which.side)]), 
          organic.GWC=mean(fungi.data$organic.GWC[which(fungi.data$Side==which.side)]),
          mineral.GWC=mean(fungi.data$mineral.GWC[which(fungi.data$Side==which.side)]),
          organic.soil.C.N=mean(fungi.data$organic.soil.C.N[which(fungi.data$Side==which.side)]),
          mineral.soil.C.N=mean(fungi.data$mineral.soil.C.N[which(fungi.data$Side==which.side)]))
fungi.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
fungi.pred1 <- as.data.frame(fungi.pred1)
fungi.pred1$Distance <- pred_interval

fungi.left <- ggplot(data=fungi.pred1, aes(x=Distance, y=fit))+ geom_line(color="red") + 
  geom_ribbon(data=fungi.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = fungi.data[which(fungi.data$Side==which.side),], aes(x = Distance, 
  y = C.N, group=Distance, color = NULL), alpha = 0.1, color="red")+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = 
  "black", fill = NA, linewidth = 1.2)) + labs(y = "C:N") + 
  theme_classic() + theme(legend.position="none", axis.title.x = element_blank(), 
  axis.text.x = element_text(size = 12),axis.title.y = element_text(size = 14),axis.text.y = element_text(size = 12)) +
  ylim(5, 14) + 
  #scale_x_continuous(breaks=c(1,3,6,10,20, 40, 60, 100)) 
  scale_x_reverse(breaks=c(1,3,6,10,20, 40, 60, 100)) 

pred_interval <- seq(1,20,0.2)
which.side <- "SR"
example <- data.frame(Distance=pred_interval, Side=which.side, 
                      carcass.deposition = "yes", 
                      total.N = mean(fungi.data$total.N[which(fungi.data$Side==which.side)]), 
                      organic.GWC=mean(fungi.data$organic.GWC[which(fungi.data$Side==which.side)]),
                      mineral.GWC=mean(fungi.data$mineral.GWC[which(fungi.data$Side==which.side)]),
                      organic.soil.C.N=mean(fungi.data$organic.soil.C.N[which(fungi.data$Side==which.side)]),
                      mineral.soil.C.N=mean(fungi.data$mineral.soil.C.N[which(fungi.data$Side==which.side)]))
fungi.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
fungi.pred1 <- as.data.frame(fungi.pred1)
fungi.pred1$Distance <- pred_interval

fungi.right <- ggplot(data=fungi.pred1, aes(x=Distance, y=fit))+ geom_line(color="red") + 
  geom_ribbon(data=fungi.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = fungi.data[which(fungi.data$Side==which.side),], aes(x = Distance, 
  y = C.N, group=Distance, color = NULL), alpha = 0.1, color="red")+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = 
  "black", fill = NA, linewidth = 1.2)) + labs(y = "C:N") + 
  theme_classic() +
  theme(panel.background = element_blank(),
        axis.text.y = element_blank(),
        axis.text.x = element_text(size = 12),
        axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        legend.position="none") +
  ylim(5, 14) +
  scale_x_continuous(breaks=c(1,3,6,10,20), expand = c(0.05, 0.05)) + 
  easy_remove_axes("y")

grid.arrange(fungi.left, fungi.right, nrow=1, bottom = textGrob("Distance",gp=gpar(fontsize=15,font=1)))

ggsave(plot=p1, width=5, height=3, dpi=300, filename = "C.N_100_fungi_L.jpg")
ggsave(plot=p2, width=5, height=3, dpi=300, filename = "C.N_100_fungi_R.jpg")

# for fungal C:N 1-100, use ylim(4.4,14)
# for fungal C:N 1-20, use ylim(5,14)

##############################################
# dN15 #
#############################################
# select data
fungi.data <- read.csv("WSU mushroom isotopes for analysis.csv")
fungi.data <- read.csv("WSU mushroom isotopes for analysis_updated_4.6.2025.csv")

# make distance numeric
fungi.data$Distance <- as.numeric(fungi.data$Distance)

# by stream 
#fungi.data <- fungi.data[ which(fungi.data$Stream=="Hansen"),]

# select only 1-20 m  
#fungi.data <- fungi.data[ which(fungi.data$Distance < 11),]

# choose trophic status
fungi.data <- fungi.data[which(fungi.data$trophic.status == "mycorrhizal" | fungi.data$trophic.status=="unknown"),]
fungi.data <- fungi.data[which(fungi.data$trophic.status == "saprotroph" | fungi.data$trophic.status=="unknown"),]

# choose exploration type
fungi.data <- fungi.data[which(fungi.data$exploration.type == "long distance" | fungi.data$trophic.status=="unknown"),]
fungi.data <- fungi.data[which(fungi.data$exploration.type == "long distance"),]
fungi.data <- fungi.data[which(fungi.data$exploration.type == "medium distance smooth"),]
fungi.data <- fungi.data[which(fungi.data$exploration.type == "medium distance fringe"),]
fungi.data <- fungi.data[which(fungi.data$exploration.type == "short distance"),]

# scale variables

test1 <- scale(fungi.data$organic.GWC, scale=FALSE)
fungi.data$organic.GWC <- test1[,1]

test2 <- scale(fungi.data$mineral.GWC, scale=FALSE)
fungi.data$mineral.GWC <- test2[,1]

test3 <- scale(fungi.data$organic.soil.dN15, scale=FALSE)
fungi.data$organic.soil.dN15 <- test3[,1]

test4 <- scale(fungi.data$mineral.soil.dN15, scale=FALSE)
fungi.data$mineral.soil.dN15 <- test4[,1]

test5 <- scale(fungi.data$organic.soil.C.N, scale=FALSE)
fungi.data$organic.soil.C.N <- test5[,1]

test6 <- scale(fungi.data$mineral.soil.C.N, scale=FALSE)
fungi.data$mineral.soil.C.N <- test6[,1]

test7 <- scale(fungi.data$X.N, scale=FALSE)
fungi.data$X.N <- test7[,1]

#test8 <- scale(fungi.data$carcass.deposition.ord, scale=FALSE)
#fungi.data$carcass.deposition.ord <- test8[,1]

fungi.data$decomposing.carcass <- as.factor(fungi.data$decomposing.carcass)
fungi.data$carcass.deposition <- as.factor(fungi.data$carcass.deposition)
fungi.data$trophic.status <- as.factor(fungi.data$trophic.status)

fungi.data$carcass.removal <- as.factor(fungi.data$carcass.removal)
fungi.data$carcass.load <- as.numeric(fungi.data$carcass.load)
fungi.data$carcass.deposition.ord <- as.factor(fungi.data$carcass.deposition.ord)

#fungi.data <- fungi.data[which(fungi.data$Distance < 11),]
#fungi.data <- fungi.data[which(fungi.data$Genus=="Cortinarius"),]

#test <- lm(d15N ~ carcass.deposition + decomposing.carcass, data=fungi.data_species)
#summary(test)

#test5 <- scale(fungi.data$C.N, scale=FALSE)
#fungi.data$C.N <- test5[,1]

# examine effect of distance on dN15 ratio

# for 1-20 with hump, omit transect and organic soil dN15:Side. For much bigger hump,
# remove fourth power interaction term

# remove certain variables (lose 13 data points for fungi) - decided not to do this because that is a lot of missing data
#fungi.data <- fungi.data[which(fungi.data$trophic.status != ""),]

# if using full dataset, remove Aleknagik data
fungi.data <- fungi.data[which(fungi.data$Stream != "Aleknagik"),]

# remove NA values for organic and mineral soil dN15 and C:N
library(tidyr)
fungi.data <- fungi.data[!is.na(fungi.data$organic.soil.dN15),] 
fungi.data <- fungi.data[!is.na(fungi.data$mineral.soil.dN15),] 
fungi.data <- fungi.data[!is.na(fungi.data$organic.soil.C.N),] 
fungi.data <- fungi.data[!is.na(fungi.data$mineral.soil.C.N),] 

# if using foilar N:P, remove any rows that do not have value for foliar N:P
#fungi.data <- fungi.data[!is.na(fungi.data$paper.birch.N.P),]

# don't lose any data for Lactarius

# change na. action
options(na.action = "na.fail")
all.parms<-lm(d15N ~ 
                Side + 
                ln(Distance) +
                Side:ln(Distance) +
                Side:I(ln(Distance)^2) + 
                #Side:I(ln(Distance)^3) + 
                #Side:I(ln(Distance)^4) + 
                X.N + 
                decomposing.carcass +
                carcass.load +
                carcass.deposition + 
                #white.spruce.N.P + 
                #paper.birch.N.P + 
                trophic.status +
                #carcass.deposition.ord + 
                #carcass.removal + 
                #ln(C.N) +
                organic.GWC + 
                organic.soil.dN15 + 
                organic.soil.C.N + 
                mineral.GWC +
                mineral.soil.dN15 +
                mineral.soil.C.N,
                #mineral.soil.dN15:Side + 
                #organic.soil.dN15:Side,
                #Genus,  
                #(1|Transect),
                data=fungi.data)
#dredge
results<-dredge(all.parms)
subset(results, delta <2)
subset(results, delta == 0)

# test model
mod1 <- lm(formula = d15N ~ ln(Distance) + X.N + decomposing.carcass + 
             carcass.load + carcass.deposition + trophic.status + organic.GWC + 
             organic.soil.dN15 + organic.soil.C.N + mineral.GWC + mineral.soil.dN15 + 
             mineral.soil.C.N, data = fungi.data)
summary(mod1)
+(1|Genus) + (1|Stream:Transect)

MuMIn::r.squaredGLMM(lmer(formula = d15N ~ ln(Distance) + X.N + decomposing.carcass + 
                            carcass.deposition + trophic.status + 
                            organic.soil.dN15 + mineral.soil.C.N + (1|Genus) + (1|Stream:Transect), data = fungi.data))

test.mod <-lmer(formula = d15N ~ ln(Distance) + X.N + decomposing.carcass + 
       carcass.deposition + trophic.status + 
       organic.soil.dN15 + mineral.soil.C.N + (1|Genus) + (1|Stream:Transect), data = fungi.data)

# long distance is distance and carcass deposition
# medium distance fringe is decomposing carcass and mineral soil d15N
# medium distance smooth is carcass deposition, distance, mineral soil C:N, organic GWC, organic soil d15N, %N
# short distance is carcass load, decomposing.carcass + ln(Distance) + mineral.GWC + organic.GWC + organic.soil.dN15 
# carcass load and mineral GWC not significant with random effects!

# step model
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

mod1 <- lmer(d15N ~ ln(Distance) + X.N + decomposing.carcass + carcass.deposition + 
               organic.GWC + organic.soil.dN15 + organic.soil.C.N + 
               trophic.status + (1|Stream:Transect) + (1|Genus), data=fungi.data)

mod2 <- lmer(d15N ~ ln(Distance) + X.N + decomposing.carcass + carcass.deposition + 
               mineral.GWC + mineral.soil.dN15 + mineral.soil.C.N + 
               trophic.status + (1|Stream:Transect) + (1|Genus), data=fungi.data)
AIC(mod1, mod2)

# with all variables
dredge.model <- lmer(d15N ~ ln(Distance) + X.N + decomposing.carcass +carcass.load + carcass.deposition + organic.GWC + 
                organic.soil.dN15 + organic.soil.C.N + mineral.GWC + mineral.soil.dN15 + mineral.soil.C.N + trophic.status +
                (1|Stream:Transect) + (1|Genus), data=fungi.data)
summary(dredge.model)

plot(dredge.model)
qqnorm(residuals(dredge.model))

# BY SPECIES
# overall
dredge.model <- lmer(d15N ~ ln(Distance) + X.N + decomposing.carcass + carcass.deposition + organic.GWC + 
                       organic.soil.dN15 + organic.soil.C.N + mineral.GWC + mineral.soil.dN15 + mineral.soil.C.N + 
                       (1|Stream:Transect), data=fungi.data)
summary(dredge.model)

# LACTARIUS
mod1 <- lmer(d15N ~ ln(Distance) + X.N + decomposing.carcass + carcass.deposition + 
               organic.GWC + organic.soil.dN15 + organic.soil.C.N + 
               (1|Stream:Transect), data=fungi.data)

mod2 <- lmer(d15N ~ ln(Distance) + X.N + decomposing.carcass + carcass.deposition + 
               mineral.GWC + mineral.soil.dN15 + mineral.soil.C.N + 
               (1|Stream:Transect), data=fungi.data)
AIC(mod1, mod2)

dredge.model <- lmer(d15N ~ carcass.deposition + ln(Distance) + 
                     organic.soil.dN15 + mineral.soil.C.N + organic.GWC + (1|Stream:Transect), data=fungi.data)
summary(dredge.model)

step.model <- lmer(d15N ~ ln(Distance) + carcass.deposition + organic.GWC + 
                     organic.soil.dN15 + mineral.soil.C.N + (1|Stream:Transect), data=fungi.data)
summary(step.model)

MuMIn::r.squaredGLMM(lme4::lmer(d15N ~ ln(Distance) + carcass.deposition + organic.GWC + 
                                  organic.soil.dN15 + mineral.soil.C.N + (1|Stream:Transect), data=fungi.data))

# all significant above

# PAXILLUS
dredge.model <- lm(d15N ~ carcass.deposition + ln(Distance) + 
                     mineral.soil.dN15 + X.N, data=fungi.data)
summary(dredge.model)

step.model <- lm(formula = d15N ~ ln(Distance) + X.N + carcass.deposition + 
    mineral.GWC + mineral.soil.dN15, data = fungi.data)

summary(step.model)


# carcass deposition

# HEBELOMA
dredge.model <- lmer(d15N ~ decomposing.carcass + organic.soil.dN15 + X.N + (1|Stream:Transect), data=fungi.data)
summary(dredge.model)

mod1 <- lm(formula = d15N ~ X.N + decomposing.carcass + organic.soil.dN15, data = fungi.data)
summary(mod1)

# LACCARIA
dredge.model <- lm(d15N ~ decomposing.carcass + carcass.deposition +
                     mineral.soil.dN15 + mineral.GWC + mineral.soil.C.N + 
                     #organic.soil.dN15 + organic.soil.C.N + organic.GWC +
                     X.N + ln(Distance), data=fungi.data)
summary(dredge.model)

mod1 <- lm(d15N ~ decomposing.carcass + carcass.deposition + carcass.load + X.N + mineral.soil.dN15, data=fungi.data)
mod2 <- lm(d15N ~ decomposing.carcass + carcass.deposition + carcass.load + X.N + organic.soil.dN15, data=fungi.data)
summary(mod2)
summary(mod1)
AIC(mod1, mod2)

# carcass deposition

# Russula
dredge.model <- lm(d15N ~ ln(Distance) + organic.soil.dN15, data=fungi.data)
summary(dredge.model)

step.model <- lmer(d15N ~ ln(Distance) + organic.GWC + organic.soil.dN15 + 
                     (1|Stream:Transect), data = fungi.data)
summary(step.model)


MuMIn::r.squaredGLMM(lme4::lmer(d15N ~ ln(Distance) + organic.GWC + organic.soil.dN15 + 
                                  (1|Stream:Transect), data = fungi.data))
# organic soil N15

# Cortinarius
dredge.model <- lm(d15N ~ ln(Distance) + organic.soil.dN15, data=fungi.data)
summary(dredge.model)

step.model <- lm(formula = d15N ~ decomposing.carcass + organic.soil.dN15 + 
     organic.soil.C.N + mineral.GWC + mineral.soil.C.N, data = fungi.data)
summary(step.model)

mod1 <- lmer(d15N ~ ln(Distance) + X.N + decomposing.carcass + 
               organic.GWC + organic.soil.dN15 + organic.soil.C.N + 
               (1|Stream:Transect), data=fungi.data)

mod2 <- lm(d15N ~ ln(Distance) + X.N + decomposing.carcass + 
               mineral.GWC + mineral.soil.dN15 + mineral.soil.C.N, data=fungi.data)

mod2 <- lm(d15N ~ ln(Distance) + X.N + decomposing.carcass + 
             mineral.GWC +  mineral.soil.C.N, data=fungi.data)


AIC(mod1, mod2)

# for all fungal genera for 1-20
dredge.model <- lmer(d15N ~ carcass.deposition + carcass.load + ln(Distance) + mineral.soil.C.N + organic.soil.dN15 + 
                        trophic.status + X.N + (1|Genus) + (1|Stream:Transect), data=fungi.data)

# for 1-20, no decomposing but organic N15

dredge.model2 <- lmer(d15N ~ carcass.deposition + decomposing.carcass + carcass.load + ln(Distance) + X.N + trophic.status + 
                      mineral.soil.C.N + mineral.soil.dN15 + (1|Genus) + (1|Stream:Transect), data=fungi.data)

# for 1-20 no decomposing but mineral C:N (this one was best)

# for all fungal genera for 1-100
dredge.model <- lmer(d15N ~ 1+ carcass.deposition + decomposing.carcass + ln(Distance) + mineral.soil.C.N + organic.soil.dN15 + 
                     trophic.status + X.N + (1|Genus) + (1|Stream:Transect), data=fungi.data)
summary(dredge.model)

# get conditional r^2 (which takes into account random effects, second number)
MuMIn::r.squaredGLMM(lme4::lmer(d15N ~ ln(Distance) + X.N + decomposing.carcass +carcass.load + carcass.deposition + trophic.status + organic.GWC + 
                                  organic.soil.dN15 + organic.soil.C.N + mineral.GWC + mineral.soil.dN15 + mineral.soil.C.N + (1|Genus) + 
                                  (1|Stream:Transect), data=fungi.data))

pred_r_squared(dredge.model)

# get percentage variance explained by each variable using anova
af <- anova(dredge.model)
afss <- af$"Sum Sq"
blah <- cbind(af,PctExp=afss/sum(afss)*100)
blah$PctExp

# calculate VIF
library(car)
car::vif(dredge.model)

# fungi 1-20
dredge.model <- lmer(d15N ~ carcass.deposition + decomposing.carcass + organic.soil.dN15 + organic.GWC + trophic.status + X.N + (1|Transect) + (1|Genus), data=fungi.data)
summary(dredge.model)

# with 1-6 m with only data with trophic status data
dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass +  trophic.status + 
                     mineral.soil.dN15, data=fungi.data)
summary(dredge.model)

# with 1-6 m with all data
dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass +trophic.status + 
                    organic.soil.dN15 + organic.GWC + X.N, data=fungi.data)
summary(dredge.model)

# for fungi 1-20 (with deposition)
dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass + organic.soil.dN15:Side + trophic.status + Side + Side:ln(Distance) + ln(Distance) + 
                     mineral.soil.dN15, data=fungi.data)
summary(dredge.model)

# for fungi 1-100 (with deposition)
dredge.model <- lmer(d15N ~ carcass.deposition + decomposing.carcass + organic.GWC + organic.soil.C.N + 
                       organic.soil.dN15 + trophic.status + X.N + (1|Genus), data=fungi.data)
summary(dredge.model)

dredge.model <- lm(d15N ~ carcass.deposition + decomposing.carcass + organic.GWC + organic.soil.dN15 + trophic.status + , data=fungi.data)
summary(dredge.model)

# for model with no side or  distance variables, 1-100
dredge.model <- lm(d15N ~ carcass.deposition + trophic.status + decomposing.carcass + organic.GWC + organic.soil.dN15:Side + X.N, 
                   data=fungi.data)
summary(dredge.model)

# for model using side and distance, 1-20
dredge.model <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + decomposing.carcass + X.N + trophic.status + organic.soil.dN15:Side, data=fungi.data)
summary(dredge.model)

# for 1-100 test with fourth power
dredge.model <- lm(d15N ~ decomposing.carcass + organic.soil.dN15 + 
                     Side + Side:I(ln(Distance)^3) + Side:I(ln(Distance)^4) + 
                     organic.soil.dN15:Side, data=fungi.data)
summary(dredge.model)

# for 1-100 with just squared (wasn't close to top model)
dredge.model <- lm(d15N ~ decomposing.carcass + ln(Distance) + organic.soil.dN15 + 
                     Side + Side:I(ln(Distance)^2) + 
                     organic.soil.dN15:Side, data=fungi.data)

# for 1-100 with just linear
dredge.model <- lm(d15N ~ decomposing.carcass + ln(Distance) + mineral.GWC + 
                     organic.soil.dN15 + Side + organic.soil.dN15:ln(Distance) + 
                     ln(Distance):Side + organic.soil.dN15:Side, data=fungi.data)

# for 1-100 dredge
dredge.model <- lm(d15N ~ decomposing.carcass + ln(Distance) + organic.soil.dN15 + Side + Side:I(ln(Distance)^2) + organic.soil.dN15:Side, data=fungi.data)
summary(dredge.model)

# for 1-20
dredge.model <- lm(d15N ~ decomposing.carcass + organic.soil.dN15 + organic.soil.dN15:Side + Side + Side:I(ln(Distance)^3), data=fungi.data)
summary(dredge.model)

step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

# calculate predictive R^2 to see if model is overfitting
pred_r_squared(step.model)

mod1 <- lm(d15N ~ Side, data=fungi.data)
mod2 <- lm(d15N ~ Side + organic.GWC, data=fungi.data)
mod2a <- lm(d15N ~ Side + mineral.GWC, data=fungi.data)
mod3 <- lm(d15N ~ ln(Distance) + organic.GWC, data=fungi.data)
mod3a <- lm(d15N ~ ln(Distance) + mineral.GWC, data=fungi.data)
mod4 <- lm(d15N ~ ln(Distance), data=fungi.data)
mod5 <- lm(d15N ~ ln(Distance) + Side, data=fungi.data)
mod6 <- lm(d15N ~ ln(Distance) + Side + organic.GWC, data=fungi.data)
mod6a <- lm(d15N ~ ln(Distance) + Side + mineral.GWC, data=fungi.data)
mod7 <- lm(d15N ~ organic.GWC, data=fungi.data)
mod8 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + organic.GWC, data=fungi.data)
mod8a <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + mineral.GWC, data=fungi.data)
mod8b <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + organic.GWC + ln(C.N), data=fungi.data)
mod9 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N, data=fungi.data)
mod10 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance), data=fungi.data)
mod11 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + organic.GWC, data=fungi.data)
mod11a <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + mineral.GWC, data=fungi.data)
mod12 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2), data=fungi.data)
mod13 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N + organic.GWC, data=fungi.data)
mod13a <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N + mineral.GWC, data=fungi.data)
mod14 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3), data=fungi.data)
mod14a <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.GWC, data=fungi.data)
mod14b <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.GWC + X.N, data=fungi.data)
mod14c <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + mineral.GWC + X.N, data=fungi.data)
mod14d <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.soil.dN15, data=fungi.data)
mod14e <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.soil.dN15 + organic.GWC, data=fungi.data)
mod14f <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.soil.dN15 + organic.GWC + ln(C.N), data=fungi.data)
mod14g <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + Side:I(ln(Distance)^4) + organic.GWC, data=fungi.data)

mod15 <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^3), data=fungi.data)
mod15a <- lm(d15N ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + Side:I(ln(Distance)^4) + organic.soil.dN15, data=fungi.data)
# although model 15s have as low AIC as model 12, they do not include individual
# variables of side and distance which are important 
mod16a <- lm(d15N ~ Side:I(ln(Distance)^2) + X.N, data=fungi.data)
mod16 <- lm(d15N ~ Side:I(ln(Distance)^3) + X.N, data=fungi.data)
mod17 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE), data=fungi.data)
mod18 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + X.N, data=fungi.data)
mod19 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + X.N + organic.GWC, data=fungi.data)
mod20 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC + organic.soil.dN15 + ln(C.N), data=fungi.data)
mod21 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + mineral.soil.dN15, data=fungi.data)
mod22 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC, data=fungi.data)
mod23 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + ln(C.N), data=fungi.data)
mod24 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC + ln(C.N), data=fungi.data)
mod25 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.soil.dN15 + ln(C.N), data=fungi.data)
mod26 <- lm(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.soil.dN15 + ln(C.N) + organic.GWC + X.N, data=fungi.data)
mod27 <- lm(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE), data=fungi.data)
mod28 <- lm(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE) + X.N, data=fungi.data)
mod29 <- lm(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE) + X.N + organic.GWC, data=fungi.data)
mod30 <- lm(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.GWC + organic.soil.dN15 + ln(C.N), data=fungi.data)
mod31 <- lm(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE) + mineral.soil.dN15, data=fungi.data)
mod32 <- lm(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.GWC, data=fungi.data)
mod33 <- lm(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE) + ln(C.N), data=fungi.data)
mod34 <- lm(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.GWC + ln(C.N), data=fungi.data)
mod35 <- lm(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.soil.dN15 + ln(C.N), data=fungi.data)
mod36 <- lm(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.soil.dN15 + ln(C.N) + organic.GWC + X.N, data=fungi.data)
mod37 <- lm(d15N ~ Side:poly(ln(Distance), 5, raw=TRUE), data=fungi.data)
mod38 <- lm(d15N ~ Side:poly(ln(Distance), 5, raw=TRUE) + organic.GWC, data=fungi.data)
mod39 <- lm(d15N ~ Side:poly(ln(Distance), 5, raw=TRUE) + organic.GWC + organic.soil.dN15 + ln(C.N), data=fungi.data)
mod40 <- lm(d15N ~ Side:poly(ln(Distance), 6, raw=TRUE), data=fungi.data)
mod41 <- lm(d15N ~ Side:poly(ln(Distance), 6, raw=TRUE) + organic.GWC, data=fungi.data)
mod42 <- lm(d15N ~ Side:poly(ln(Distance), 6, raw=TRUE) + organic.GWC + organic.soil.dN15 + ln(C.N), data=fungi.data)

mod17.lmer <- lmer(d15N ~ Side:poly(ln(Distance), 3, raw=TRUE) + (1|Transect), data=fungi.data, REML=FALSE)
mod35.lmer <- lmer(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.soil.dN15 + ln(C.N) + (1|Transect), data=fungi.data, REML=FALSE)
mod32.lmer <- lmer(d15N ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.GWC + (1|Transect), data=fungi.data)

mod32 <- lm(d15N ~ poly(ln(Distance), 4, raw=TRUE) + organic.GWC, data=fungi.data)

results <- AIC(mod1, mod2, mod2a, mod3, mod3a, mod4, mod5, mod6, mod6a, mod7, mod8, mod8a, 
    mod8b,mod9, mod10, mod11, mod11a, mod12, mod13, mod13a, mod14, mod14a,
    mod14b,mod14c, mod14d, mod14e, mod14f, mod14g, mod15, mod15a, mod16, mod16a, mod17, mod18, mod19, 
    mod20, mod21, mod22, mod23, mod24, mod25, mod26, mod27, mod28, mod29, mod30,
    mod31, mod32, mod33, mod34, mod35, mod36, mod37, mod38, mod39, mod40, mod41,
    mod42, step.model, dredge.model)
# results for one side only
#results <- AIC(mod3, mod3a, mod4, mod7, mod32, dredge.model, step.model)
ordered <- results[order(-results$AIC),]
ordered
summary(test.model)

# FUNGI 1-100
# chose dredge


# step model shows hump on both sides, while dredge model doesn't show hump
# bigger left hump with manual model selection shows left hump for fourth best model 
# (mod32)...so could probably go either way. Mod32 doesn't have any significant 
# variables though. When running sides separately, for left and right the step shows 
# hump, but dredge doesn't show hump for either side. These are top models. 

# 

# for 1-100, many models in 30s and 20s work well (similar AIC), choosing 
# between 35, 33, 28, 25, 23/21/17/14/14d, which means that the covariates 
# include total N, C:N, organic soil dN15 and organic GWC; I chose mod32

# FUNGI 1-20
# chose step

step.model <- lm(d15N ~ Side + ln(Distance) + X.N + decomposing.carcass + organic.soil.dN15 + Side:I(ln(Distance)^3) + Side:I(ln(Distance)^4) + 
                          Side:organic.soil.dN15, data=fungi.data)

step.model.lmer <- lmer(d15N ~ Side + ln(Distance) + X.N + decomposing.carcass + 
    organic.soil.dN15 + Side:I(ln(Distance)^3) + Side:I(ln(Distance)^4) + 
      Side:organic.soil.dN15 + (1|Genus), data=fungi.data)

# select the model that will be used to plot predictions
select.mod <- dredge.model

pred_interval <- c("yes", "no")
example <- data.frame(carcass.deposition=pred_interval, 
                      X.N = mean(fungi.data$X.N), 
                      organic.GWC=mean(fungi.data$organic.GWC),
                      organic.soil.C.N=mean(fungi.data$organic.soil.C.N),
                      organic.soil.dN15=mean(fungi.data$organic.soil.dN15),
                      decomposing.carcass="no",
                      trophic.status = "mycorrhizal",
                      Side="SL")
fungi.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
fungi.pred1 <- as.data.frame(fungi.pred1)
fungi.pred1$carcass.deposition <- pred_interval

ggplot(data=fungi.pred1, aes(x=carcass.deposition, y=fit))+ 
  geom_point(color = "red") +
  geom_errorbar(data=fungi.pred1, aes(ymin=lwr, ymax=upr),width = 0.2) +
  geom_point(data = fungi.data, aes(x = carcass.deposition, y = d15N, color="red"), alpha = 0.3)+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = "black", fill = NA, linewidth = 1.2)) + 
  labs(y = "d15N", x="Carcass Deposition") + theme_classic()+theme(legend.position="none")

# plot raw data with model estimates for LEFT bank
pred_interval <- seq(1,20,0.2)
#pred_interval <- c(1,6,10,20)
which.side <- "SL"
example <- data.frame(Distance=pred_interval, Side=which.side, 
           X.N = mean(fungi.data$X.N[which(fungi.data$Side==which.side)]), 
           organic.GWC=mean(fungi.data$organic.GWC[which(fungi.data$Side==which.side)]),
           organic.soil.C.N=mean(fungi.data$organic.soil.C.N[which(fungi.data$Side==which.side)]),
  organic.soil.dN15=mean(fungi.data$organic.soil.dN15[which(fungi.data$Side==which.side)]),
  mineral.soil.dN15=mean(fungi.data$mineral.soil.dN15[which(fungi.data$Side==which.side)]),
  decomposing.carcass="yes",
  carcass.deposition="no", 
  trophic.status = "saprotroph")
fungi.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
fungi.pred1 <- as.data.frame(fungi.pred1)
fungi.pred1$Distance <- pred_interval

p1 <- ggplot(data=fungi.pred1, aes(x=Distance, y=fit))+ 
   geom_line(color = "red") +
  #geom_errorbar(data=fungi.pred1, aes(ymin=lwr, ymax=upr),width = 0.2) +
  geom_ribbon(data=fungi.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = fungi.data[which(fungi.data$Side==which.side),], aes(x = Distance, 
  y = d15N, group=Distance, color="red"), alpha = 0.3)+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = 
   "black", fill = NA, linewidth = 1.2)) + 
  labs(y = "d15N", x="Distance") + theme_classic()+theme(legend.position="none")+
  ylim(-4.88, 23.5) + 
  scale_x_reverse(breaks=c(1,3,6,10,20)) 

# plot raw data with model estimates for RIGHT bank
pred_interval <- seq(1,20,0.2)
#pred_interval <- c(1,6,10,20)
which.side <- "SR"
example <- data.frame(Distance=pred_interval, Side=which.side, 
                      X.N = mean(fungi.data$X.N[which(fungi.data$Side==which.side)]), 
                      organic.GWC=mean(fungi.data$organic.GWC[which(fungi.data$Side==which.side)]),
                      mineral.GWC=mean(fungi.data$mineral.GWC[which(fungi.data$Side==which.side)]),
                      C.N=mean(fungi.data$C.N[which(fungi.data$Side==which.side)]),
                      organic.soil.dN15=mean(fungi.data$organic.soil.dN15[which(fungi.data$Side==which.side)]),
                      mineral.soil.dN15=mean(fungi.data$mineral.soil.dN15[which(fungi.data$Side==which.side)]),
                      organic.soil.C.N=mean(fungi.data$organic.soil.C.N[which(fungi.data$Side==which.side)]),
                      decomposing.carcass="yes",
                      carcass.deposition="no",
                      trophic.status = "saprotroph")
fungi.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
fungi.pred1 <- as.data.frame(fungi.pred1)
fungi.pred1$Distance <- pred_interval

p2 <- ggplot(data=fungi.pred1, aes(x=Distance, y=fit))+ 
  geom_line(color = "red") +
  #geom_errorbar(data=fungi.pred1, aes(ymin=lwr, ymax=upr),width = 0.2) +
  geom_ribbon(data=fungi.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = fungi.data[which(fungi.data$Side==which.side),], aes(x = Distance, 
  y = d15N, group=Distance, color="red"), alpha = 0.3)+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = 
  "black", fill = NA, linewidth = 1.2)) + 
  labs(y = "d15N", x="Distance") + theme_classic()+theme(legend.position="none")+
  ylim(-4.88, 23.5) + 
  scale_x_continuous(breaks=c(1,3,6,10,20, 40, 60, 100)) 

grid.arrange(p1, p2, nrow=1)

ggsave(plot=p1, width=5, height=3, dpi=300, filename = "dN15_100_fungi_left.jpg")
ggsave(plot=p2, width=5, height=3, dpi=300, filename = "dN15_100_fungi_right.jpg")

# for 1-20 m data, use ylim(1.3, 22)
# for 1-100 m data, use ylim(-4.88, 23.5)

#############################################
# dC13 #
#############################################
# select data
fungi.data <- read.csv("WSU mushroom isotopes for analysis.csv")

# make distance numeric
fungi.data$Distance <- as.numeric(fungi.data$Distance)

# by stream 
fungi.data <- fungi.data[ which(fungi.data$Stream=="Hansen"),]

# select only 1-20 m  
fungi.data <- fungi.data[ which(fungi.data$Distance < 21),]

# scale variables
test <- scale(fungi.data$X.N, scale=FALSE)
fungi.data$X.N <- test[,1]

test1 <- scale(fungi.data$organic.GWC, scale=FALSE)
fungi.data$organic.GWC <- test1[,1]

test2 <- scale(fungi.data$mineral.GWC, scale=FALSE)
fungi.data$mineral.GWC <- test2[,1]

test3 <- scale(fungi.data$organic.soil.dC13, scale=FALSE)
fungi.data$organic.soil.dC13 <- test3[,1]

test4 <- scale(fungi.data$mineral.soil.dC13, scale=FALSE)
fungi.data$mineral.soil.dC13 <- test4[,1]

test5 <- scale(fungi.data$white.spruce.dC13, scale=FALSE)
fungi.data$white.spruce.dC13 <- test5[,1]

test6 <- scale(fungi.data$paper.birch.dC13, scale=FALSE)
fungi.data$paper.birch.dC13 <- test6[,1]

# examine effect of distance on dC13 ratio
# change na. action

fungi.data$decomposing.carcass <- as.factor(fungi.data$decomposing.carcass)
fungi.data$carcass.deposition <- as.factor(fungi.data$carcass.deposition)
fungi.data$trophic.status <- as.factor(fungi.data$trophic.status)
fungi.data$carcass.load <- as.numeric(fungi.data$carcass.load)

# if using full dataset, remove Aleknagik data
fungi.data <- fungi.data[which(fungi.data$Stream != "Aleknagik"),]

# remove certain variables (lose 13 data points for fungi)
#fungi.data <- fungi.data[which(fungi.data$trophic.status == "mycorrhizal"),]

# include mineral GWC and mineral soil dC13? Not sure...
library(tidyr)

# if including plant dC13, need to drop a few rows that did not have data for plant dC13 
fungi.data <- fungi.data %>% drop_na(white.spruce.dC13)
fungi.data <- fungi.data %>% drop_na(paper.birch.dC13)

options(na.action = "na.fail")
all.parms<-lm(dC13 ~ 
                #Side + 
                ln(Distance) + 
                #Side:ln(Distance) +
                #Side:I((ln(Distance))^2) + 
                #Side:I((ln(Distance))^3) + 
                #Side:I((ln(Distance))^4) + 
                #ln(C.N) + 
                white.spruce.dC13 + 
                #Side:white.spruce.dC13 +
                paper.birch.dC13 +
                #Side:paper.birch.dC13 +
                X.N + 
                organic.soil.dC13 + 
                mineral.GWC + 
                mineral.soil.dC13 +
                decomposing.carcass + 
                carcass.deposition + 
                trophic.status + 
                organic.GWC,
              data=fungi.data)
results<-dredge(all.parms)
subset(results, delta <2)
subset(results, delta == 0)

# all data
dredge.model <- lmer(dC13 ~ ln(Distance) + white.spruce.dC13 + paper.birch.dC13 + X.N + organic.soil.dC13 + mineral.soil.dC13 +
                     decomposing.carcass + carcass.deposition + trophic.status + (1|Genus) + (1|Stream:Transect), data=fungi.data)
summary(dredge.model)

dredge.model <- lm(dC13 ~ ln(Distance) + mineral.soil.dC13 + organic.GWC + paper.birch.dC13 + trophic.status + X.N, data=fungi.data)
summary(dredge.model)

dredge.model <- lmer(dC13 ~ ln(Distance) + mineral.soil.dC13 + organic.GWC + paper.birch.dC13 + trophic.status + X.N + 
                       (1|Stream:Transect), data=fungi.data)
summary(dredge.model)

# get conditional r^2 (which takes into account random effects, second number)
MuMIn::r.squaredGLMM(lme4::lmer(dC13 ~ ln(Distance) + mineral.soil.dC13 + organic.GWC + paper.birch.dC13 + trophic.status + X.N + 
                                  (1|Stream:Transect), data=fungi.data))

# get percentage variance explained by each variable using anova
af <- anova(dredge.model)
afss <- af$"Sum Sq"
blah <- cbind(af,PctExp=afss/sum(afss)*100)
blah$PctExp

# calculate VIF
library(car)
car::vif(dredge.model)

# for 1-20 fungi dC13 (with deposition variable)
dredge.model <- lm(dC13 ~ paper.birch.dC13 + ln(Distance) + X.N, data=fungi.data)
summary(dredge.model)

# for 1-100 fungi dC13 (with deposition variable)
dredge.model <- lm(dC13 ~ X.N + trophic.status, data=fungi.data)
summary(dredge.model)

# for 1-100
dredge.model <- lm(dC13 ~ ln(Distance) + X.N, data=fungi.data)
summary(dredge.model)

# for 1-20
dredge.model <- lm(dC13 ~ ln(Distance) + X.N + paper.birch.dC13, data=fungi.data)
summary(dredge.model)

# stepwise regression
step.model <- stepAIC(all.parms, direction = "both", trace = FALSE)
summary(step.model)

pred_r_squared(mod26)

mod1 <- lm(dC13 ~ Side, data=fungi.data)
mod2 <- lm(dC13 ~ Side + organic.GWC, data=fungi.data)
mod2a <- lm(dC13 ~ Side + mineral.GWC, data=fungi.data)
mod3 <- lm(dC13 ~ ln(Distance) + organic.GWC, data=fungi.data)
mod3a <- lm(dC13 ~ ln(Distance) + mineral.GWC, data=fungi.data)
mod4 <- lm(dC13 ~ ln(Distance), data=fungi.data)
mod5 <- lm(dC13 ~ ln(Distance) + Side, data=fungi.data)
mod6 <- lm(dC13 ~ ln(Distance) + Side + organic.GWC, data=fungi.data)
mod6a <- lm(dC13 ~ ln(Distance) + Side + mineral.GWC, data=fungi.data)
mod7 <- lm(dC13 ~ organic.GWC, data=fungi.data)
mod8 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + organic.GWC, data=fungi.data)
mod8a <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + mineral.GWC, data=fungi.data)
mod8b <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2) + organic.GWC + ln(C.N), data=fungi.data)
mod9 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N, data=fungi.data)
mod10 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance), data=fungi.data)
mod11 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + organic.GWC, data=fungi.data)
mod11a <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + mineral.GWC, data=fungi.data)
mod12 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I((ln(Distance))^2), data=fungi.data)
mod13 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N + organic.GWC, data=fungi.data)
mod13a <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + X.N + mineral.GWC, data=fungi.data)
mod14 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3), data=fungi.data)
mod14a <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.GWC, data=fungi.data)
mod14b <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.GWC + X.N, data=fungi.data)
mod14c <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + mineral.GWC + X.N, data=fungi.data)
mod14d <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.soil.dC13, data=fungi.data)
mod14e <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.soil.dC13 + organic.GWC, data=fungi.data)
mod14f <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^2) + Side:I(ln(Distance)^3) + organic.soil.dC13 + organic.GWC + ln(C.N), data=fungi.data)

mod15 <- lm(dC13 ~ Side + ln(Distance) + Side:ln(Distance) + Side:I(ln(Distance)^3), data=fungi.data)
# although model 15s have as low AIC as model 12, they do not include individual
# variables of side and distance which are important 
mod16a <- lm(dC13 ~ Side:I(ln(Distance)^2) + X.N, data=fungi.data)
mod16 <- lm(dC13 ~ Side:I(ln(Distance)^3) + X.N, data=fungi.data)
mod17 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE), data=fungi.data)
mod18 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + X.N, data=fungi.data)
mod19 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + X.N + organic.GWC, data=fungi.data)
mod20 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC + organic.soil.dC13 + ln(C.N), data=fungi.data)
mod21 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + mineral.soil.dC13, data=fungi.data)
mod22 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC, data=fungi.data)
mod23 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + ln(C.N), data=fungi.data)
mod24 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.GWC + ln(C.N), data=fungi.data)
mod25 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.soil.dC13 + ln(C.N), data=fungi.data)
mod26 <- lm(dC13 ~ Side:poly(ln(Distance), 3, raw=TRUE) + organic.soil.dC13 + ln(C.N) + organic.GWC + X.N, data=fungi.data)
mod27 <- lm(dC13 ~ Side:poly(ln(Distance), 4, raw=TRUE), data=fungi.data)
mod28 <- lm(dC13 ~ Side:poly(ln(Distance), 4, raw=TRUE) + X.N, data=fungi.data)
mod29 <- lm(dC13 ~ Side:poly(ln(Distance), 4, raw=TRUE) + X.N + organic.GWC, data=fungi.data)
mod30 <- lm(dC13 ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.GWC + organic.soil.dC13 + ln(C.N), data=fungi.data)
mod31 <- lm(dC13 ~ Side:poly(ln(Distance), 4, raw=TRUE) + mineral.soil.dC13, data=fungi.data)
mod32 <- lm(dC13 ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.GWC, data=fungi.data)
mod33 <- lm(dC13 ~ Side:poly(ln(Distance), 4, raw=TRUE) + ln(C.N), data=fungi.data)
mod34 <- lm(dC13 ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.GWC + ln(C.N), data=fungi.data)
mod35 <- lm(dC13 ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.soil.dC13 + ln(C.N), data=fungi.data)
mod36 <- lm(dC13 ~ Side:poly(ln(Distance), 4, raw=TRUE) + organic.soil.dC13 + ln(C.N) + organic.GWC + X.N, data=fungi.data)
mod37 <- lm(dC13 ~ Side:poly(ln(Distance), 5, raw=TRUE), data=fungi.data)
mod38 <- lm(dC13 ~ Side:poly(ln(Distance), 5, raw=TRUE) + organic.GWC, data=fungi.data)
mod39 <- lm(dC13 ~ Side:poly(ln(Distance), 5, raw=TRUE) + organic.GWC + organic.soil.dC13 + ln(C.N), data=fungi.data)
mod40 <- lm(dC13 ~ Side:poly(ln(Distance), 6, raw=TRUE), data=fungi.data)
mod41 <- lm(dC13 ~ Side:poly(ln(Distance), 6, raw=TRUE) + organic.GWC, data=fungi.data)
mod42 <- lm(dC13 ~ Side:poly(ln(Distance), 6, raw=TRUE) + organic.GWC + organic.soil.dC13 + ln(C.N), data=fungi.data)

results <- AIC(mod1, mod2, mod2a, mod3, mod3a, mod4, mod5, mod6, mod6a, mod7, mod8, mod8a, 
    mod8b,mod9, mod10, mod11, mod11a, mod12, mod13, mod13a, mod14, mod14a, 
    mod14b,mod14c, mod14d, mod14e, mod14f, mod15, mod16, mod16a, mod17, mod18, mod19, 
    mod20, mod21, mod22, mod23, mod24, mod25, mod26, mod27, mod28, mod29, mod30,
    mod31, mod32, mod33, mod34, mod35, mod36, mod37, mod38, mod39, mod40, mod41,
    mod42, step.model, dredge.model)
ordered <- results[order(-results$AIC),]
ordered
summary(mod14b)

# for fungal dC13 1-100, best model is dredge
# for fungal dC13 1-20, chose step

# for fungal dC13 1-100, use ylim(-29.6, -21.8)
# for fungal dC13 1-20, use  ylim(-30, -23.5) 

# select the model that will be used to plot predictions
select.mod <- step.model

# plot raw data with model estimates
pred_interval <- seq(1,20,0.2)
which.side <- "SL"
example <- data.frame(Distance=pred_interval, Side=which.side, 
           X.N = mean(fungi.data$X.N[which(fungi.data$Side==which.side)]), 
           organic.GWC=mean(fungi.data$organic.GWC[which(fungi.data$Side==which.side)]),
           mineral.GWC=mean(fungi.data$mineral.GWC[which(fungi.data$Side==which.side)]),
           organic.soil.dC13=mean(fungi.data$organic.soil.dC13[which(fungi.data$Side==which.side)]),
           mineral.soil.dC13=mean(fungi.data$mineral.soil.dC13[which(fungi.data$Side==which.side)]),
           C.N=mean(fungi.data$C.N[which(fungi.data$Side==which.side)]),
           paper.birch.dC13=mean(fungi.data$paper.birch.dC13[which(fungi.data$Side==which.side)]))
fungi.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
fungi.pred1 <- as.data.frame(fungi.pred1)
fungi.pred1$Distance <- pred_interval

p1 <- ggplot(data=fungi.pred1, aes(x=Distance, y=fit))+ geom_line(color="red") + 
  geom_ribbon(data=fungi.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = fungi.data[which(fungi.data$Side==which.side),], aes(x = Distance, 
  y = dC13, group=Distance, color = NULL), alpha = 0.1, color="red")+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = 
  "black", fill = NA, linewidth = 1.2)) + labs(y = "d13C", x="Distance") + 
  theme_classic()+ theme(legend.position="none") +
  ylim(-30, -23.5) + 
  scale_x_reverse(breaks=c(1,3,6,10,20, 40, 60, 100)) 

# plot raw data with model estimates for RIGHT BANK
pred_interval <- seq(1,20,0.2)
which.side <- "SR"
example <- data.frame(Distance=pred_interval, Side=which.side, 
                      X.N = mean(fungi.data$X.N[which(fungi.data$Side==which.side)]), 
                      organic.GWC=mean(fungi.data$organic.GWC[which(fungi.data$Side==which.side)]),
                      mineral.GWC=mean(fungi.data$mineral.GWC[which(fungi.data$Side==which.side)]),
                      organic.soil.dC13=mean(fungi.data$organic.soil.dC13[which(fungi.data$Side==which.side)]),
                      mineral.soil.dC13=mean(fungi.data$mineral.soil.dC13[which(fungi.data$Side==which.side)]),
                      C.N=mean(fungi.data$C.N[which(fungi.data$Side==which.side)]),
                      paper.birch.dC13=mean(fungi.data$paper.birch.dC13[which(fungi.data$Side==which.side)]))
fungi.pred1 <- predict(select.mod, newdata=example, re.form=NA, type="response", interval='confidence')
fungi.pred1 <- as.data.frame(fungi.pred1)
fungi.pred1$Distance <- pred_interval

p2 <- ggplot(data=fungi.pred1, aes(x=Distance, y=fit))+ geom_line(color="red") + 
  geom_ribbon(data=fungi.pred1, aes(ymin=lwr, ymax=upr), linetype=0, alpha=0.1) + 
  geom_point(data = fungi.data[which(fungi.data$Side==which.side),], aes(x = Distance, 
  y = dC13, group=Distance, color = NULL), alpha = 0.1, color="red")+
  theme(panel.background = element_blank(),panel.border = element_rect(colour = 
  "black", fill = NA, linewidth = 1.2)) + labs(y = "d13C", x="Distance") + 
  theme_classic() + theme(legend.position="none") +
  ylim(-30, -23.5)  + 
  scale_x_continuous(breaks=c(1,3,6,10,20, 40, 60, 100)) 

grid.arrange(p1, p2, nrow=1)

ggsave(plot=p1, width=5, height=3, dpi=300, filename = "dC13_20_fungi_L.jpg")
ggsave(plot=p2, width=5, height=3, dpi=300, filename = "dC13_20_fungi_R.jpg")

##############################################
# %N #
#########################################
# select data
fungi.data <- read.csv("WSU mushroom isotopes for analysis.csv")

# make distance numeric
fungi.data$Distance <- as.numeric(fungi.data$Distance)

# by stream 
#fungi.data <- fungi.data[ which(fungi.data$Stream=="Hansen"),]

# select only 1-20 m  
#fungi.data <- fungi.data[ which(fungi.data$Distance < 21),]

# scale variables
test1 <- scale(fungi.data$organic.GWC, scale=FALSE)
fungi.data$organic.GWC <- test1[,1]

test2 <- scale(fungi.data$mineral.GWC, scale=FALSE)
fungi.data$mineral.GWC <- test2[,1]

test3 <- scale(fungi.data$organic.soil.C.N, scale=FALSE)
fungi.data$organic.soil.C.N <- test3[,1]

test4 <- scale(fungi.data$mineral.soil.C.N, scale=FALSE)
fungi.data$mineral.soil.C.N <- test4[,1]

fungi.data$decomposing.carcass <- as.factor(fungi.data$decomposing.carcass)
fungi.data$carcass.deposition <- as.factor(fungi.data$carcass.deposition)
fungi.data$trophic.status <- as.factor(fungi.data$trophic.status)
fungi.data$carcass.deposition.ord <- as.factor(fungi.data$carcass.deposition.ord)
fungi.data$carcass.load <- as.numeric(fungi.data$carcass.load)

# if using full dataset, remove Aleknagik data
fungi.data <- fungi.data[which(fungi.data$Stream != "Aleknagik"),]

# remove NA values for organic and mineral soil dN15 and C:N
library(tidyr)
fungi.data <- fungi.data[!is.na(fungi.data$organic.soil.C.N),] 
fungi.data <- fungi.data[!is.na(fungi.data$mineral.soil.C.N),] 

# remove certain variables (lose 13 data points for fungi)
#fungi.data <- fungi.data[which(fungi.data$trophic.status == "mycorrhizal"),]

# if including plant C:N, need to drop a few rows that did not have data for plant C:N 
fungi.data <- fungi.data %>% drop_na(white.spruce.C.N)
fungi.data <- fungi.data %>% drop_na(paper.birch.C.N)

fungi.data <- fungi.data[which(fungi.data$Genus=="Laccaria"),]

# parameters for BOTH SIDES
options(na.action = "na.fail")
all.parms<-lm(X.N ~ 
                #Side + 
                ln(Distance) + 
                #Side:ln(Distance) +
                #Side:I(ln(Distance)^2) + 
                #Side:I(ln(Distance)^3) + 
                #Side:I(ln(Distance)^4) + 
                paper.birch.C.N + 
                white.spruce.C.N +
                organic.GWC + 
                mineral.GWC +
                organic.soil.C.N + 
                mineral.soil.C.N + 
                trophic.status + 
                carcass.deposition +
                #carcass.load + 
                #carcass.deposition.ord + 
                decomposing.carcass,
              data=fungi.data)
results<-dredge(all.parms)
subset(results, delta <2)
subset(results, delta == 0)

# all data
dredge.model <- lmer(X.N ~ ln(Distance) + paper.birch.C.N + white.spruce.C.N + organic.GWC + mineral.GWC + organic.soil.C.N + 
                       mineral.soil.C.N + trophic.status + carcass.deposition + decomposing.carcass + carcass.load + (1|Genus) + 
                       (1|Stream:Transect), data=fungi.data)
summary(dredge.model)

# from dredge
dredge.model <- lmer(X.N ~ carcass.deposition + paper.birch.C.N + 
                       trophic.status + white.spruce.C.N + (1|Genus) + 
                       (1|Stream:Transect), data=fungi.data)

# model for single genera
dredge.model <- lm(C.N ~ carcass.deposition + decomposing.carcass, data=fungi.data)
summary(dredge.model)

dredge.model <- lm(C.N ~ carcass.deposition.ord + paper.birch.C.N + trophic.status + white.spruce.C.N, data=fungi.data)
summary(dredge.model)

# model for all generas from 1-100
dredge.model <- lm(C.N ~ carcass.deposition.ord + ln(Distance) + paper.birch.C.N + trophic.status + white.spruce.C.N, data=fungi.data)
summary(dredge.model)

dredge.model <- lmer(C.N ~ carcass.deposition + ln(Distance) + paper.birch.C.N + trophic.status + white.spruce.C.N + 
                       (1|Genus), data=fungi.data)
summary(dredge.model)

# get conditional r^2 (which takes into account random effects, second number)
MuMIn::r.squaredGLMM(lme4::lmer(C.N ~ ln(Distance) + paper.birch.C.N + white.spruce.C.N + organic.GWC + mineral.GWC + organic.soil.C.N + 
                                  mineral.soil.C.N + trophic.status + carcass.deposition + decomposing.carcass + carcass.load + (1|Genus) + 
                                  (1|Stream:Transect), data=fungi.data))
